library(tidyverse)
library(magrittr)
library(parallel)
library(ngsReports)
library(here)
library(scales)
library(ggpubr)
library(kableExtra)
library(AnnotationHub)
library(ensembldb)
library(edgeR)
library(corrplot)
library(DT)
library(ggrepel)
library(msigdbr)
library(fgsea)
if (interactive()) setwd(here::here())
theme_set(theme_bw())
cores <- detectCores() - 2
ah_Dr <- AnnotationHub() %>%
subset(species == "Danio rerio") %>%
subset(rdataclass == "EnsDb")
ensDb <- ah_Dr[["AH83189"]]
trEns_Dr <- transcripts(ensDb) %>%
mcols() %>%
as_tibble()
trLen_Dr <- exonsBy(ensDb, "tx") %>%
width() %>%
vapply(sum, integer(1))
geneGcLen_Dr <- trLen_Dr %>%
enframe() %>%
set_colnames(c("tx_id", "length")) %>%
left_join(trEns_Dr) %>%
group_by(gene_id) %>%
summarise(
aveLen = mean(length),
maxLen = max(length),
aveGc = sum(length * gc_content) / sum(length),
longestGc = gc_content[which.max(length)[[1]]]
) %>%
mutate(
aveGc = aveGc / 100,
longestGc = longestGc / 100
)
trGcLen_Dr <- trLen_Dr %>%
enframe() %>%
set_colnames(c("tx_id", "length")) %>%
left_join(trEns_Dr) %>%
group_by(tx_id) %>%
summarise(
aveLen = mean(length),
maxLen = max(length),
aveGc = sum(length * gc_content) / sum(length),
longestGc = gc_content[which.max(length)[[1]]]
) %>%
mutate(
aveGc = aveGc / 100,
longestGc = longestGc / 100
)
genesGR_Dr <- genes(ensDb)
mcols(genesGR_Dr) <- mcols(genesGR_Dr)[c("gene_id", "gene_name",
"gene_biotype", "entrezid")]
txGR_Dr <- transcripts(ensDb)
mcols(txGR_Dr) <- mcols(txGR_Dr)[c("tx_id", "tx_name",
"tx_biotype", "tx_id_version", "gene_id")]
An EnsDb object was obtained for Ensembl release 101
using the AnnotationHub package. This provided the GC
content and length for every gene and transcript in the release. For
zebrafish, this consists of 37241 genes and 65905 transcripts.
This is a total RNA-seq dataset generated from a 3-way comparison of WT zebrafish (Danio rerio) with heterozygous mutants (psen2S4Ter/+) and homozygous mutants (psen2S4Ter/S4Ter). A previous analysis of this dataset identified the possibility of incomplete ribosomal RNA (rRNA) removal. The following analysis involves an investigation into possible reasons for incomplete rRNA removal and any bias this introduces into the data.
files <- list.files(
path = "/mnt/hpcfs/a1647910/20200310_rRNADepletion/1_Psen2S4Ter/0_rawData/FastQC",
pattern = "zip",
full.names = TRUE
)
samples <- tibble(
sample = str_remove(basename(files), "_fastqc.zip"),
dataset = NA,
organism = NA
) %>%
mutate(
dataset = ifelse(
str_detect(sample, "Ps2Ex"), "Psen2S4Ter", dataset
),
organism = ifelse(
str_detect(sample, "Ps2Ex"), "zebrafish", organism
)
)
datasets <- samples$dataset %>%
unique()
The following analysis involves 12 paired-end samples across 1 dataset(s): Psen2S4Ter.
rawFqc <- files %>%
FastqcDataList()
data <- grep("Ps2Ex", fqName(rawFqc))
labels <- rawFqc[data] %>%
fqName() %>%
str_remove("_6month_07_07_2016_F3") %>%
str_remove("\\.fastq\\.gz") %>%
str_remove("Ps2Ex3M1_")
rawLib <- plotReadTotals(rawFqc[data]) +
labs(subtitle = "Psen2S4Ter") +
scale_x_discrete(labels = labels)
The library sizes of the unprocessed dataset(s) range between 27,979,654 and 37,144,975 reads.
rawLib
rRNA transcripts are known to have high GC content. Therefore, inspecting the GC content of the raw reads is a logical start point for detecting incomplete rRNA removal. A spike in GC content at > 70% is expected if this is the case.
plotly::ggplotly(
plotGcContent(
x = rawFqc[data],
plotType = "line",
gcType = "Transcriptome",
species = "Drerio"
) +
labs(title = "Psen2S4Ter Dataset (D. rerio)") +
theme(legend.position="none")
)
GC content of reads in the dataset. Clear spikes above 70% GC are observed, which is likely due to incomplete rRNA depletion.
The top 30 overrepresented sequences were analysed using
blastn and were found to be predominantly rRNA
sequences.
getModule(rawFqc, "Overrep") %>%
group_by(Sequence, Possible_Source) %>%
summarise(`Found In` = n(), `Highest Percentage` = max(Percentage)) %>%
arrange(desc(`Highest Percentage`), desc(`Found In`)) %>%
ungroup() %>%
dplyr::slice(1:30) %>%
mutate(`Highest Percentage` = percent_format(0.01)(`Highest Percentage`/100)) %>%
kable(
align = "llrr",
caption = paste(
"Top", nrow(.),"Overrepresented sequences.",
"The number of samples they were found in is shown,",
"along with the percentage of the most 'contaminated' sample."
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive")
)
| Sequence | Possible_Source | Found In | Highest Percentage |
|---|---|---|---|
| GTGGGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGC | No Hit | 12 | 1.72% |
| CTGGGGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAG | No Hit | 10 | 1.69% |
| CCGCTGTATTACTCAGGCTGCACTGCAGTGTCTATTCACAGGCGCGATCC | No Hit | 12 | 1.32% |
| GGCCCGGCGCACGTCCAGAGTCGCCGCCGCACACCGCAGCGCATCCCCCC | No Hit | 9 | 1.31% |
| CTCCTGAAAAGGTTGTATCCTTTGTTAAAGGGGCTGTACCCTCTTTAACT | No Hit | 11 | 1.11% |
| GGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCATC | No Hit | 12 | 1.09% |
| GGGGTGTACGAAGCTGAACTTTTATTCATCTCCCAGACAACCAGCTATTG | No Hit | 12 | 1.07% |
| GGCCCGGCGCACGTCCAGAGTCGCCGCCGCGCACCGCAGCGCATCCCCCC | No Hit | 10 | 1.03% |
| CGAGAGGCTCTAGTTGATATACTACGGCGTAAAGGGTGGTTAAGGAACAA | No Hit | 12 | 1.01% |
| GGGGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAGCG | No Hit | 9 | 0.87% |
| CCTCCTTCAAGTATTGTTTCATGTTACATTTTCGTATATTCTGGGGTAGA | No Hit | 12 | 0.82% |
| CCCGCTGTATTACTCAGGCTGCACTGCAGTGTCTATTCACAGGCGCGATC | No Hit | 12 | 0.79% |
| GTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCATCT | No Hit | 12 | 0.79% |
| GGGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCAT | No Hit | 12 | 0.78% |
| CGGGTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGG | No Hit | 12 | 0.73% |
| GGGCCTCTAGCATCTAAAAGCGTATAACAGTTAAAGGGCCGTTTGGCTTT | No Hit | 11 | 0.68% |
| CAGTGGCGTGCGCCTGTAATCCAAGCTACTGGGAGGCTGAGGCTGGCGGA | No Hit | 11 | 0.64% |
| CGGGTCGGGTGGGTAGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGG | No Hit | 12 | 0.57% |
| CTTAGACGACCTGGTAGTCCAAGGCTCCCCCAGGAGCACCATATCGATAC | No Hit | 11 | 0.54% |
| AGCTGGGGAGATCCGCGAGAAGGGCCCGGCGCACGTCCAGAGTCGCCGCC | No Hit | 11 | 0.53% |
| GGCCTCTAGCATCTAAAAGCGTATAACAGTTAAAGGGCCGTTTGGCTTTA | No Hit | 10 | 0.53% |
| CAGCCTATTTAACTTAGGGCCAACCCGTCTCTGTGGCAATAGAGTGGGAA | No Hit | 12 | 0.51% |
| GGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGACGTGG | No Hit | 10 | 0.50% |
| GGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAGCGCG | No Hit | 9 | 0.49% |
| CTGGGAGATGAATAAAAGTTCAGCTTCGTACACCCCAAATTAAAAAATTA | No Hit | 10 | 0.48% |
| GCCTATTTAACTTAGGGCCAACCCGTCTCTGTGGCAATAGAGTGGGAAGA | No Hit | 12 | 0.47% |
| GGTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGAC | No Hit | 11 | 0.45% |
| CCCCCGAACCCTTCCAAGCCGAACCGGAGCCGGTCGCGGCGCACCGCCGA | No Hit | 10 | 0.45% |
| GTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGACG | No Hit | 10 | 0.43% |
| GCCCACTACGACAACGTGTTTTGTAAATTATGATCTTTATTCTCCTGAAA | No Hit | 10 | 0.43% |
Raw libraries were trimmed using cutadapt v1.14 to
remove Illumina adapter sequences. Bases with PHRED score < 30,
NextSeq-induced polyG runs and reads shorter than 35bp were also
removed.
trimFqc <- list.files(
path = "/mnt/hpcfs/a1647910/20200310_rRNADepletion/1_Psen2S4Ter/1_trimmedData/FastQC",
pattern = "zip",
full.names = TRUE
) %>%
FastqcDataList()
trimStats <- readTotals(rawFqc) %>%
dplyr::rename(Raw = Total_Sequences) %>%
left_join(readTotals(trimFqc), by = "Filename") %>%
dplyr::rename(Trimmed = Total_Sequences) %>%
mutate(
Discarded = 1 - Trimmed/Raw,
Retained = Trimmed / Raw
)
After trimming of adapters between 4.16% and 5.08% of reads were discarded.
Trimmed reads were:
Aligned to rRNA sequences using the BWA-MEM
algorithm to estimate the proportion of reads that were of rRNA origin
within each sample. BWA-MEM is recommended for high-quality
queries of reads ranging from 70bp to 1Mbp as it is faster and more
accurate that alternative algorithms BWA-backtrack and
BWA-SW.
Aligned to the Danio rerio GRCz11 genome (Ensembl
release 101) using STAR v2.7.0d and summarised with
featureCounts from the Subread v1.5.2 package.
These counts were used for all gene-level analysis.
rRnaProp <- read.delim(
"/mnt/hpcfs/a1647910/20200310_rRNADepletion/1_Psen2S4Ter/3_bwa/log/samples.mapped.all",
sep = ":",
col.names = c("sample", "proportion"),
header = FALSE
) %>%
mutate(
sample = str_remove_all(sample, "_6month_F3|[0-9]*_Ps2Ex3M1_|.mapped"),
sample = basename(sample),
proportion = proportion/100,
dataset = "Psen2S4Ter",
organism = "zebrafish",
group = str_extract(sample, "(WT|Heter|Hom)")
) %>%
as_tibble()
rRnaProp$dataset %<>%
factor(levels = c("Psen2S4Ter"))
rRnaProp %>%
ggplot(aes(sample, proportion)) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~dataset, scales = "free_x") +
scale_y_continuous(labels = percent) +
labs(x = "Sample", y = "Percent of Total", fill = "Read pair") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Percentages of each library that align to rRNA sequences with
bwa mem.
rRnaProp %>%
ggplot(aes(group, proportion, fill = group)) +
geom_boxplot() +
scale_y_continuous(labels = percent) +
labs(x = "Genotype", y = "Percent of total RNA", title = "rRNA proportions of each genotype") +
scale_fill_discrete(
name = "Genotype"
)
dgeList <- read_tsv("/mnt/hpcfs/a1647910/20200310_rRNADepletion/1_Psen2S4Ter/4_star2pass/featureCounts/genes.out") %>%
set_colnames(basename(colnames(.))) %>%
set_colnames(str_remove(colnames(.), "Aligned.sortedByCoord.out.bam")) %>%
set_colnames(str_remove(colnames(.), "_6month_F3")) %>%
set_colnames(str_remove(colnames(.), "[0-9]*_Ps2Ex3M1_")) %>%
as.data.frame() %>%
column_to_rownames("Geneid") %>%
DGEList() %>%
calcNormFactors()
dgeList$genes <- genesGR_Dr[rownames(dgeList),]
mcols(dgeList$genes) %<>%
as.data.frame() %>%
left_join(geneGcLen_Dr)
addInfo <- tibble(
sample = rRnaProp$sample,
dataset = "Psen2S4Ter",
organism = "zebrafish",
rRNA = rRnaProp$proportion
)
dgeList$samples %<>%
rownames_to_column("rowname") %>%
mutate(sample = rowname) %>%
left_join(addInfo) %>%
column_to_rownames("rowname")
# dgeList$samples$filenames <- list.files(
# "/mnt/hpcfs/a1647910/20200310_rRNADepletion/1_Psen2S4Ter/2_alignedData/bam",
# pattern = ".bam$",
# full.names = TRUE
# )
dgeList$samples$group <- colnames(dgeList) %>%
str_extract("(WT|Heter|Hom)") %>%
factor(levels = c("WT", "Heter", "Hom"))
gcInfo <- function(x) {
x$counts %>%
as.data.frame() %>%
rownames_to_column("gene_id") %>%
as_tibble() %>%
pivot_longer(
cols = colnames(x),
names_to = "sample",
values_to = "counts"
) %>%
dplyr::filter(
counts > 0
) %>%
left_join(
geneGcLen_Dr
) %>%
dplyr::select(
ends_with("id"), sample, counts, aveGc, maxLen
) %>%
split(f = .$sample) %>%
lapply(
function(x){
DataFrame(
gc = Rle(x$aveGc, x$counts),
logLen = Rle(log10(x$maxLen), x$counts)
)
}
)
}
gcSummary <- function(x) {
x %>%
vapply(function(x){
c(mean(x$gc), sd(x$gc), mean(x$logLen), sd(x$logLen))
}, numeric(4)
) %>%
t() %>%
set_colnames(
c("mn_gc", "sd_gc", "mn_logLen", "sd_logLen")
) %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
as_tibble()
}
rle <- gcInfo(dgeList)
sumGc <- gcSummary(rle)
a <- sumGc %>%
left_join(dgeList$samples) %>%
ggplot(aes(rRNA, mn_logLen)) +
geom_point(aes(colour = group), size = 3) +
geom_smooth(method = "lm") +
scale_x_continuous(labels = percent) +
labs(
x = "rRNA Proportion of Initial Library",
y = "Mean log(Length)",
colour = "Genotype"
)
b <- sumGc %>%
left_join(dgeList$samples) %>%
ggplot(aes(rRNA, mn_gc)) +
geom_point(aes(colour = group), size = 3) +
geom_smooth(method = "lm") +
scale_y_continuous(labels = percent) +
scale_x_continuous(labels = percent) +
labs(
x = "rRNA Proportion of Initial Library",
y = "Mean GC Content",
colour = "Genotype"
)
ggarrange(
a, b, ncol = 2, nrow = 1,
common.legend = TRUE, legend = "bottom"
) %>%
annotate_figure("PsenS4Ter Dataset (D. rerio)")
Comparison of residual bias potentially introduced by incomplete rRNA removal. Regression lines are shown along with standard error bands for each comparison.
genes2keep <- dgeList %>%
cpm() %>%
is_greater_than(1) %>%
rowSums() %>%
is_weakly_greater_than(6)
dgeFilt <- dgeList[genes2keep,, keep.lib.sizes = FALSE] %>%
calcNormFactors()
pca <- cpm(dgeFilt, log = TRUE) %>%
t() %>%
prcomp()
pcaCor <- pca$x %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sumGc) %>%
as_tibble() %>%
left_join(dgeList$samples) %>%
dplyr::select(
PC1, PC2, PC3,
Mean_GC = mn_gc,
Mean_Length = mn_logLen,
rRna_Proportion = rRNA
) %>%
cor()
a <- pca$x %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(dgeList$samples) %>%
as_tibble() %>%
ggplot(aes(PC1, PC2)) +
geom_point(aes(colour = group), size = 2) +
labs(
x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
y = paste0("PC2 (", percent(summary(pca)$importance["Proportion of Variance","PC2"]),")"),
colour = "Genotype"
)
b <- pca$x %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(dgeList$samples) %>%
ggplot(aes(PC1, rRNA, label = rRNA)) +
geom_point(aes(colour = group), size = 2) +
geom_smooth(method = "lm") +
geom_text_repel(show.legend = FALSE) +
scale_y_continuous(labels = percent) +
labs(
x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
y = "rRNA Proportion of Initial Library",
colour = "Genotype"
)
c <- pca$x %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sumGc) %>%
left_join(dgeList$samples) %>%
as_tibble() %>%
ggplot(aes(PC1, mn_logLen)) +
geom_point(aes(colour = group), size = 2) +
geom_smooth(method = "lm") +
labs(
x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
y = "Mean log(Length)",
colour = "Genotype"
)
d <- pca$x %>%
as.data.frame() %>%
rownames_to_column("sample") %>%
left_join(sumGc) %>%
left_join(dgeList$samples) %>%
as_tibble() %>%
ggplot(aes(PC1, mn_gc)) +
geom_point(aes(colour = group), size = 2) +
geom_smooth(method = "lm") +
scale_y_continuous(labels = percent) +
labs(
x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
y = "Mean GC",
colour = "Genotype"
)
ggarrange(
a, b, c, d, ncol = 2, nrow = 2,
common.legend = TRUE, legend = "bottom"
) %>%
annotate_figure("Psen2S4Ter")
PCA plot showing rRNA proportion, mean GC content and mean log(length) after summarisation to gene-level.
corrplot(
pcaCor,
type = "lower",
diag = FALSE,
addCoef.col = 1, addCoefasPercent = TRUE
)
Correlations between the first three principal components and measured variables: mean GC content, mean log(length) and rRNA proportion.
displayRes_de <- function(x){
de <- x %>%
dplyr::filter(DE)
de %>%
dplyr::slice(1:1000) %>%
dplyr::select(-gene_biotype, -coef, -DE) %>%
mutate(across(c("P.Value", "FDR", "Bonf"), ~ sprintf("%.2e", .x))) %>%
datatable(
options = list(pageLength = 10),
class = "striped hover condensed responsive",
filter = "top",
caption = paste0(
x$coef[1],
": ",
nrow(de),
" of ",
nrow(x),
" genes were classified as differentially expressed ",
"with a FDR < 0.05. ",
"If more than 1000 genes were classified as DE, only the top 1000 are shown."
)
) %>%
formatRound(c("logFC", "logCPM", "F"), digits = 2)
}
design_r <- model.matrix(~rRNA, data = dgeFilt$samples)
disp_r <- estimateDisp(dgeFilt, design_r)
plotBCV(disp_r)
fit_r <- glmQLFit(disp_r)
plotQLDisp(fit_r)
res_r <- glmQLFTest(fit_r) %>%
topTags(n = Inf) %>%
.[["table"]] %>%
rename_all(
str_remove, pattern = "ID."
) %>%
dplyr::select(
Geneid = gene_id, Symbol = gene_name, gene_biotype, logFC, logCPM, F,
P.Value = PValue, FDR, aveLen, maxLen, aveGc, longestGc
) %>%
as_tibble() %>%
mutate(
Bonf = p.adjust(P.Value, "bonf"),
DE = FDR < 0.05
)
res_r %>%
dplyr::filter(DE) %>%
dplyr::slice(1:1000) %>%
dplyr::select(-gene_biotype, -DE, -aveLen, -maxLen, -aveGc, -longestGc) %>%
mutate(across(c("P.Value", "FDR", "Bonf"), ~ sprintf("%.2e", .x))) %>%
datatable(
options = list(pageLength = 10),
class = "striped hover condensed responsive",
filter = "top",
caption = paste0(
"rRNA: ",
nrow(dplyr::filter(res_r, DE)),
" of ",
nrow(res_r),
" genes were classified as differentially expressed ",
"with a FDR < 0.05. ",
"If more than 1000 genes were classified as DE, only the top 1000 are shown."
)
) %>%
formatRound(c("logFC", "logCPM", "F"), digits = 2)
design_g <- model.matrix(~group, data = dgeFilt$samples)
disp_g <- estimateDisp(dgeFilt, design_g)
plotBCV(disp_g)
fit_g <- glmQLFit(disp_g)
plotQLDisp(fit_g)
res_g <- colnames(design_g) %>%
str_subset("Intercept", negate = TRUE) %>%
sapply(function(x){
glmQLFTest(fit_g, coef = x) %>%
topTags(n = Inf) %>%
.[["table"]] %>%
rename_all(
str_remove, pattern = "ID."
) %>%
dplyr::select(
Geneid = gene_id, Symbol = gene_name, gene_biotype, logFC, logCPM, F,
P.Value = PValue, FDR
) %>%
as_tibble() %>%
mutate(
Bonf = p.adjust(P.Value, "bonf"),
coef = x,
DE = FDR < 0.05
)
},
simplify = FALSE)
displayRes_de(res_g$groupHeter)
displayRes_de(res_g$groupHom)
# voom <- voom(dgeFilt, design = design)
# fit <- lmFit(voom, design = design)
# eBayes <- eBayes(fit)
# topTable <- eBayes %>%
# topTable(coef = colnames(design)[2], sort.by = "p", n = Inf) %>%
# set_colnames(str_remove(colnames(.), "ID\\.")) %>%
# mutate(Bonf = p.adjust(P.Value, "bonferroni")) %>%
# mutate(DE = Bonf < 0.05) %>%
# unite(Location, c(seqnames, start, end, width, strand), sep = ":") %>%
# dplyr::select(
# Geneid = gene_id,
# Symbol = gene_name,
# AveExpr,
# logFC,
# P.Value,
# FDR = adj.P.Val,
# Location,
# t,
# DE,
# everything(),
# -B
# ) %>%
# as_tibble()
design_gr <- model.matrix(~rRNA + group, data = dgeFilt$samples)
disp_gr <- estimateDisp(dgeFilt, design_gr)
plotBCV(disp_gr)
fit_gr <- glmQLFit(disp_gr)
plotQLDisp(fit_gr)
res_gr <- colnames(design_gr) %>%
str_subset("Intercept", negate = TRUE) %>%
sapply(function(x){
glmQLFTest(fit_gr, coef = x) %>%
topTags(n = Inf) %>%
.[["table"]] %>%
rename_all(
str_remove, pattern = "ID."
) %>%
dplyr::select(
Geneid = gene_id, Symbol = gene_name, gene_biotype, logFC, logCPM, F,
P.Value = PValue, FDR,
) %>%
as_tibble() %>%
mutate(
Bonf = p.adjust(P.Value, "bonf"),
coef = x,
DE = FDR < 0.05
)
},
simplify = FALSE)
displayRes_de(res_gr$rRNA)
displayRes_de(res_gr$groupHeter)
displayRes_de(res_gr$groupHom)
design_deconv <- design_gr %>%
as.data.frame() %>%
rownames_to_column("ID") %>%
mutate(
groupWT = as.integer(str_detect(ID, "WT")),
across(
.cols = starts_with("group"),
.fns = function(x){
x <- x - .$rRNA
x[x<0] <- 0
x
}
)
) %>%
dplyr::select(
ID, rRNA, starts_with("group")
) %>%
column_to_rownames("ID") %>%
as.matrix()
cont_deconv <-makeContrasts(
groupHeter = groupHeter - groupWT,
groupHom = groupHom - groupWT,
rRNA = rRNA,
levels = colnames(design_deconv)
)
disp_deconv <- estimateDisp(dgeFilt, design_deconv)
plotBCV(disp_deconv)
fit_deconv <- glmQLFit(disp_deconv)
plotQLDisp(fit_deconv)
res_deconv <- colnames(cont_deconv) %>%
sapply(function(x){
glmQLFTest(fit_deconv, contrast = cont_deconv[,x]) %>%
topTags(n = Inf) %>%
.[["table"]] %>%
rename_all(
str_remove, pattern = "ID."
) %>%
dplyr::select(
Geneid = gene_id, Symbol = gene_name, gene_biotype, logFC, logCPM, F,
P.Value = PValue, FDR,
) %>%
as_tibble() %>%
mutate(
Bonf = p.adjust(P.Value, "bonf"),
coef = x,
DE = FDR < 0.05
)
},
simplify = FALSE)
displayRes_de(res_deconv$rRNA)
displayRes_de(res_deconv$groupHeter)
displayRes_de(res_deconv$groupHom)
entrezGenes <- mcols(dgeList$genes) %>%
as.data.frame() %>%
dplyr::filter(!is.na(entrezid)) %>%
unnest(cols = entrezid) %>%
dplyr::rename(entrez_gene = entrezid)
ranks_r <- res_r %>%
mutate(stat = -sign(logFC) * log10(P.Value)) %>%
dplyr::arrange(stat) %>%
with(structure(stat, names = Geneid))
ranks_g <- res_g %>%
lapply(function(x){
x %>%
mutate(stat = -sign(logFC) * log10(P.Value)) %>%
dplyr::arrange(stat) %>%
with(structure(stat, names = Geneid))
})
ranks_gr <- res_gr %>%
lapply(function(x){
x %>%
mutate(stat = -sign(logFC) * log10(P.Value)) %>%
dplyr::arrange(stat) %>%
with(structure(stat, names = Geneid))
})
ranks_deconv <- res_deconv %>%
lapply(function(x){
x %>%
mutate(stat = -sign(logFC) * log10(P.Value)) %>%
dplyr::arrange(stat) %>%
with(structure(stat, names = Geneid))
})
displayRes_enrich <- function(x, cap){
x %>%
unnest(cols = leadingEdge) %>%
group_by(pathway) %>%
mutate(
leadingSize = n(),
pathway = str_remove(pathway, "HALLMARK_|KEGG_|WP_"),
pathway = str_trunc(pathway, 33)
) %>%
distinct(pathway, .keep_all = TRUE) %>%
dplyr::select(-leadingEdge) %>%
mutate(across(c("pval", "FDR", "padj"), ~ sprintf("%.2e", .x))) %>%
datatable(
options = list(pageLength = 10),
class = "striped hover condensed responsive",
filter = "top",
caption = paste(cap)
) %>%
formatRound(c("log2err", "ES", "NES"), digits = 2)
}
hm <- msigdbr("Danio rerio", category = "H") %>%
left_join(entrezGenes) %>%
dplyr::filter(!is.na(gene_id)) %>%
distinct(gs_name, gene_id, .keep_all = TRUE)
hmByGene <- hm %>%
split(f = .$gene_id) %>%
lapply(extract2, "gs_name")
hmByID <- hm %>%
split(f = .$gs_name) %>%
lapply(extract2, "gene_id")
kg <- msigdbr("Danio rerio", category = "C2", subcategory = "CP:KEGG") %>%
left_join(entrezGenes) %>%
dplyr::filter(!is.na(gene_id)) %>%
distinct(gs_name, gene_id, .keep_all = TRUE)
kgByGene <- kg %>%
split(f = .$gene_id) %>%
lapply(extract2, "gs_name")
kgByID <- kg %>%
split(f = .$gs_name) %>%
lapply(extract2, "gene_id")
wk <- msigdbr("Danio rerio", category = "C2", subcategory = "CP:WIKIPATHWAYS") %>%
left_join(entrezGenes) %>%
dplyr::filter(!is.na(gene_id)) %>%
distinct(gs_name, gene_id, .keep_all = TRUE)
wkByGene <- wk %>%
split(f = .$gene_id) %>%
lapply(extract2, "gs_name")
wkByID <- wk %>%
split(f = .$gs_name) %>%
lapply(extract2, "gene_id")
gsSizes <- bind_rows(hm, kg, wk) %>%
dplyr::select(gs_name, gene_symbol, gene_id) %>%
chop(c(gene_symbol, gene_id)) %>%
mutate(
gs_size = vapply(gene_symbol, length, integer(1))
)
fgsea_r_hm <-fgsea(hmByID, ranks_r, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
displayRes_enrich(fgsea_r_hm, "rRNA")
fgsea_g_hm <- ranks_g %>%
lapply(function(x){
fgsea(hmByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_g_hm$groupHeter, "groupHeter")
displayRes_enrich(fgsea_g_hm$groupHom, "groupHom")
fgsea_gr_hm <- ranks_gr %>%
lapply(function(x){
fgsea(hmByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_gr_hm$rRNA, "rRNA")
displayRes_enrich(fgsea_gr_hm$groupHeter, "groupHeter")
displayRes_enrich(fgsea_gr_hm$groupHom, "groupHom")
fgsea_deconv_hm <- ranks_deconv %>%
lapply(function(x){
fgsea(hmByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_deconv_hm$rRNA, "rRNA")
displayRes_enrich(fgsea_deconv_hm$groupHeter, "groupHeter")
displayRes_enrich(fgsea_deconv_hm$groupHom, "groupHom")
fgsea_r_kg <-fgsea(kgByID, ranks_r, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
displayRes_enrich(fgsea_r_kg, "rRNA")
fgsea_g_kg <- ranks_g %>%
lapply(function(x){
fgsea(kgByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_g_kg$groupHeter, "groupHeter")
displayRes_enrich(fgsea_g_kg$groupHom, "groupHom")
fgsea_gr_kg <- ranks_gr %>%
lapply(function(x){
fgsea(kgByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_gr_kg$rRNA, "rRNA")
displayRes_enrich(fgsea_gr_kg$groupHeter, "groupHeter")
displayRes_enrich(fgsea_gr_kg$groupHom, "groupHom")
fgsea_deconv_kg <- ranks_deconv %>%
lapply(function(x){
fgsea(kgByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_deconv_kg$rRNA, "rRNA")
displayRes_enrich(fgsea_deconv_kg$groupHeter, "groupHeter")
displayRes_enrich(fgsea_deconv_kg$groupHom, "groupHom")
fgsea_r_wk <-fgsea(wkByID, ranks_r, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
displayRes_enrich(fgsea_r_wk, "rRNA")
fgsea_g_wk <- ranks_g %>%
lapply(function(x){
fgsea(wkByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_g_wk$groupHeter, "groupHeter")
displayRes_enrich(fgsea_g_wk$groupHom, "groupHom")
fgsea_gr_wk <- ranks_gr %>%
lapply(function(x){
fgsea(wkByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_gr_wk$rRNA, "rRNA")
displayRes_enrich(fgsea_gr_wk$groupHeter, "groupHeter")
displayRes_enrich(fgsea_gr_wk$groupHom, "groupHom")
fgsea_deconv_wk <- ranks_deconv %>%
lapply(function(x){
fgsea(wkByID, x, eps = 0) %>%
as_tibble() %>%
dplyr::rename(FDR = padj) %>%
mutate(padj = p.adjust(pval, "bonferroni")) %>%
dplyr::arrange(pval)
})
displayRes_enrich(fgsea_deconv_wk$rRNA, "rRNA")
displayRes_enrich(fgsea_deconv_wk$groupHeter, "groupHeter")
displayRes_enrich(fgsea_deconv_wk$groupHom, "groupHom")
save(
addInfo,
dgeList,
dgeFilt,
res_r,
res_g,
res_gr,
res_deconv,
file = here::here(
"1_Psen2S4Ter/R/output/1_1_DE.RData"
)
)
sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Australia/Adelaide
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] fgsea_1.32.0 msigdbr_7.5.1 ggrepel_0.9.6
## [4] DT_0.33 corrplot_0.95 edgeR_4.4.0
## [7] limma_3.62.1 ensembldb_2.28.0 AnnotationFilter_1.28.0
## [10] GenomicFeatures_1.56.0 AnnotationDbi_1.66.0 Biobase_2.66.0
## [13] GenomicRanges_1.56.1 GenomeInfoDb_1.40.1 IRanges_2.38.1
## [16] S4Vectors_0.42.1 AnnotationHub_3.12.0 BiocFileCache_2.14.0
## [19] dbplyr_2.5.0 kableExtra_1.4.0 ggpubr_0.6.0
## [22] scales_1.3.0 here_1.0.1 ngsReports_2.7.1
## [25] patchwork_1.3.0 BiocGenerics_0.52.0 magrittr_2.0.3
## [28] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [31] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [34] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [37] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 ggdendro_0.2.0
## [3] rstudioapi_0.17.1 jsonlite_1.8.9
## [5] farver_2.1.2 rmarkdown_2.28
## [7] BiocIO_1.14.0 zlibbioc_1.52.0
## [9] vctrs_0.6.5 memoise_2.0.1
## [11] Rsamtools_2.20.0 RCurl_1.98-1.16
## [13] rstatix_0.7.2 htmltools_0.5.8.1
## [15] S4Arrays_1.4.1 curl_6.0.1
## [17] broom_1.0.7 SparseArray_1.4.8
## [19] sass_0.4.9 bslib_0.8.0
## [21] htmlwidgets_1.6.4 plotly_4.10.4
## [23] zoo_1.8-12 cachem_1.1.0
## [25] GenomicAlignments_1.40.0 mime_0.12
## [27] lifecycle_1.0.4 pkgconfig_2.0.3
## [29] Matrix_1.7-1 R6_2.5.1
## [31] fastmap_1.2.0 GenomeInfoDbData_1.2.13
## [33] MatrixGenerics_1.18.0 digest_0.6.37
## [35] colorspace_2.1-1 rprojroot_2.0.4
## [37] crosstalk_1.2.1 RSQLite_2.3.8
## [39] labeling_0.4.3 filelock_1.0.3
## [41] fansi_1.0.6 timechange_0.3.0
## [43] mgcv_1.9-1 httr_1.4.7
## [45] abind_1.4-8 compiler_4.4.2
## [47] bit64_4.5.2 withr_3.0.2
## [49] backports_1.5.0 BiocParallel_1.40.0
## [51] carData_3.0-5 DBI_1.2.3
## [53] ggsignif_0.6.4 MASS_7.3-61
## [55] rappdirs_0.3.3 DelayedArray_0.30.1
## [57] rjson_0.2.23 tools_4.4.2
## [59] glue_1.8.0 restfulr_0.0.15
## [61] nlme_3.1-166 grid_4.4.2
## [63] generics_0.1.3 gtable_0.3.6
## [65] tzdb_0.4.0 data.table_1.16.2
## [67] hms_1.1.3 xml2_1.3.6
## [69] car_3.1-2 utf8_1.2.4
## [71] XVector_0.44.0 BiocVersion_3.20.0
## [73] pillar_1.9.0 vroom_1.6.5
## [75] babelgene_22.9 splines_4.4.2
## [77] lattice_0.22-6 rtracklayer_1.64.0
## [79] bit_4.5.0 tidyselect_1.2.1
## [81] locfit_1.5-9.10 Biostrings_2.72.1
## [83] knitr_1.49 gridExtra_2.3
## [85] ProtGenerics_1.38.0 SummarizedExperiment_1.34.0
## [87] svglite_2.1.3 xfun_0.49
## [89] statmod_1.5.0 matrixStats_1.4.1
## [91] stringi_1.8.4 UCSC.utils_1.0.0
## [93] lazyeval_0.2.2 yaml_2.3.10
## [95] evaluate_1.0.1 codetools_0.2-20
## [97] BiocManager_1.30.25 cli_3.6.3
## [99] systemfonts_1.1.0 munsell_0.5.1
## [101] jquerylib_0.1.4 Rcpp_1.0.13-1
## [103] png_0.1-8 XML_3.99-0.17
## [105] blob_1.2.4 bitops_1.0-9
## [107] viridisLite_0.4.2 crayon_1.5.3
## [109] rlang_1.1.4 fastmatch_1.1-4
## [111] cowplot_1.1.3 KEGGREST_1.44.1